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1. Introduction 

The physics of disordered quantum systems is a very active field of research since the 
1950s, when this topic received the fundamental contributions of Anderson and Mott [I]. 
Their works showed that, in presence of disorder, even ideal conductors undergo a phase 
transition towards an insulating state due to destructive quantum interference. An ideal 
playground for the investigation of quantum effects is offered by ultracold atoms, where 
a controllable amount of disorder may be implemented in many ways, e.g., by means 
of speckle potentials, bichromatic lattices with incommensurate frequencies, localized 
impurities, or site-resolved addressing in optical lattices. After a long search, Anderson 
localization of ultracold ideal gases was finally observed in one dimension (ID) j2], and 
very recently reported also in 3D (3j. 

At this point it is important to stress the fundamental difference between disordered 
fermionic and bosonic systems. For fermions the crucial role is played by the Fermi 
energy: the question whether a given system is a conductor or insulator depends on 
whether the single particle states close to Fermi level are localized or not. Bosons, in 
contrast, tend to condense in the lowest energy state (or states). The question whether 
they constitute a superfluid or an insulator reduces then to understanding if the low- 
energy states are localized or not (for more extensive discussion see for instance [4j|5]). 

Dirty bosons phase diagram /su+vyt 




u/t 

Interaction strength 



Figure 1. Schematic phase diagram for dirty bosons at zero temperature as a function 
of interactions (U) and chemical potential ([i) containing glassy (G), superfluid (BEC), 
and Mott insulating (MI) phases. The continuous and dashed lines mark respectively 
the superfluid-insulator and glass-Mott transition, while the dotted line indicates a 
crossover from the Lifshits glass (heavily fragmented) to the Bose glass (extended) 
region. Outside of the BEC region, the one-body density matrix p(r,r') decays as 
|r — r'| — > oo. The Mott insulator lobes appear at sufficiently strong interactions, as 
the inclusion theorem guarantees that there may be no MI for U < V. 
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In the past decades there appeared an enormous literature on disordered systems, 
see, e.g., (4 14 , although due to the complexity of the topic mathematically rigorous 



results are still scarce. A great deal of work has been done on random Schrodinger 
operators, which describe single particle behaviour, see 15 for an excellent survey. In 
the recent paper 16 two of us have provided rigorous tight bounds on the ground state 
energy, as well as approximations of the ground state and excited states wave-functions 
for the case of random impurity (Bernoulli) disorder in a ID lattice. Since in this 
case the shape of the low-energy states is known, there are analytical ways to estimate 
interaction energy and obtain rigorous results also for interacting systems. 

A challenging and still open problem is the understanding of interaction effects on 
disordered systems. A crucial step forward in our understanding of dirty bosonic systems 
was made by Giamarchi and Schulz |6) and Fisher et al. [7j who conjectured that, in 
presence of weak disorder, interactions give rise to an intermediate compressible but 
insulating phase, the Bose glass, in between the superfluid and the insulating strongly- 
interacting phases (Tonks gas in continuous ID systems, or Mott insulator in lattice 
systems). The disputed controversy, as to whether or not the insulating Mott phase is 
always surrounded by a glassy phase was finally settled in [13] , who demonstrated, using 
the theorem of inclusions, that this is indeed the case for the dirty boson system with 
generic bounded disorder. Experimentally, this problem has been recently addressed 
in 



17-19 



In the regime where disorder dominates over interactions, it was shown instead 
that the ground state of the quantum fluid may be described in terms of a qualitatively 
different glassy phase, called Lifshits (or Anderson) glass 10 20 , characterized by 
exponentially localized and well separated "islands". The latter may be identified 
with the low-lying single-particle eigenstates residing in the regions of space where the 
potential is small. As repulsive interactions are strengthened, an increasing number of 
islands is populated, until their overlap becomes sufficient to establish phase coherence 
and transport between them. The g whole then undergoes a phase transition 

towards an extended Bose-Einstein condensate (BEC). If the gas is confined in a lattice, 
for even larger interactions the repulsion between particles becomes so strong that 
on-site density fluctuations characteristic of the superfluid state become energetically 
unfavourable. In this situation the gas undergoes a further transition out of the BEC, 
first to the Bose glass phase and then into the incompressible Mott state. Between the 
Lifshits and the Bose glasses there is no phase transition, as both are gapless, insulating 
and compressible phases, but nonetheless the two states are qualitatively different, in 
the sense that in the former the gas is fragmented into independent and distant low- 
energy islands (the Lifshits states), while the latter tends to extend over a large portion 
of the available volume. The phases discussed above are sketched in figure [TJ and their 
most important properties listed in table [TJ 

In the present paper we consider disorder of Bernoulli type, i.e., randomly- 
distributed, identical and localized impurities, and we focus on the regime of weak 
to intermediate interactions in presence of such impurities. We first investigate the 
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properties of the ideal gas, discussing the typical size and energy of the localized islands, 
and showing that the low-energy states indeed form a so-called Lifshits tail. 

We further explore interaction effects on the ground state properties by comparing 
the performance of two theoretical methods. The first one is a simple version of 
the multi-orbital Hartree-Fock method (sMOHF), a variational method based on an 
expansion on single-particle states. This approximation is appropriate to describe 



the glassy regime where superfluidity is suppressed by disorder 10 . More elaborate 
versions of the Hartree-Fock method, that incorporate self-consistency and allow to 
describe fragmentation and quantum dynamics of interacting Bose systems can be found 
in [2l] and references therein. The second method is the standard Gross-Pitaevskii 
(GP) approach, which generally yields an appropriate description for every interaction 
strength. 

We calculate the ground state energy of the system, the associated superfluid 
fraction, the fractal dimension, and introduce a new parameter, the fractional 
occupation, allowing one to discern between the Lifshits and the Bose glass. Finally, 
we turn to the investigation of temperature effects in an ideal gas. Due to the increase 
of kinetic energy, temperature yields effects similar to those generated by repulsive 
interactions, in the sense that the gas occupies a larger number of islands, which 
increasingly overlap until the ground state covers a large amount of the available space. 
Nonetheless, we will point out that there are important differences between the two 
cases. 

As an important result, in this paper we show that disordered systems with 
Bernoulli potential allow for analytical estimates which are very well supported by 
numerical simulations. These results provided us with intuitions on how to generalize 
the rigorous results of (16) to non-interacting 2D systems, and, at least partially, to 
interacting ID and 2D systems. These rigorous results go beyond the scope of the present 
article, and will be published in a more appropriate mathematical physics journal. 

This paper is structured as follows. We introduce the disordered potential under 
study and the relevant Hamiltonian in Sec. [2] In Sec. [3] we identify the eigenstates of 
the ideal system, we numerically demonstrate the Lifshits tail behaviour, and we discuss 
the ground state energy scaling as a function of the size of the system. In Sec. [4] we 
introduce the sMOHF and GP methods, and in Sec. [5] we compute a number of relevant 
quantities for interacting systems. In Sec. [6] we compare the effects of interactions and 
non-zero temperature, and we present our conclusions in Sec. [7| 



Table 1. Summary of most common phases for dirty bosons in a lattice, and associated 



properties. 
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2. The model 

We study the properties of a bosonic gas on disordered ID and 2D square lattices 
of linear dimension L. We consider disorder of the so-called Bernoulli type, i.e., the 
potential on each site is an independent random variable that assumes the value: 

I y>0 with probability 1 — p , ^ 

| with probability p 

This kind of potential may be ideally realized by using a two-component gas, and a 



species-selective lattice which deeply traps only one of the two components 19 , 22 
At sufficiently low energies, the component which is not trapped experiences s-wave 
collisions against localized scatterers. Alternatively, a Bernoulli disorder may be 



imprinted directly on the gas exploiting single-site addressable optical lattices 23 . 

The use of Bernoulli disorder is particularly appealing because its asymptotic 
properties, as discussed in the next sections, become apparent even for small lattices 



e.g., with 50 sites). As we show in Appendix A, the Bernoulli potential has actually 



optimal properties of convergence. Therefore, the Bernoulli disorder is suitable for 
numerical simulations and due to its simple form, it allows for various analytical 
estimates. 

The Hamiltonian of the interacting system we consider is then 

M i 

where t denotes the tunnelling energy, c] is the creation operator of bosons on site i, 
and denotes nearest neighbour pairs of sites. The interaction term U int will be 

discussed later in Section HI 

3. Properties of a non-interacting gas 

Let us start by discussing the properties of a non-interacting system, and in particular 
look at its ground state energy, and the density of low energy excited states. 

3.1. Ground state energy 

In a large ID system with Bernoulli disorder, the linear size Lmax of the largest island 



of contiguous zero-potential sites scales as l£2 oc logL 16 . The proof of this fact 
is based on the following simple observation. Let Iq be the characteristic size of the 
largest island, then p l ° estimates its occurrence. The number of islands of similar size 
is expected to be of order of L (strictly speaking Lj logL as we shall discuss below), so 
that the probability of occurring of any one of them is Lp l ° — > const as L — > oo, which 
indeed implies Z — logL/| \og(p)\. It is then clear that the ground state energy of the 
non-interacting bosonic gas in such a potential will be bounded above by the energy 
of the first harmonic, i.e., a half-sine function, of the largest island. A lower bound 
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of the same order was proven in [16] , showing that the ground state energy behaves 
asymptotically as 7r 2 /(logL) 2 when L — > oo. 

Similarly, one can argue that in a large D-dimensional system the ground state 
will be occupying the largest spherical island of zero-potential sites since this shape 
minimizes the kinetic energy [24]. Its diameter can be shown to grow as L^ x oc 
(\ogL) l l D , and therefore we expect the energy of the ground state to scale as E oc 
logL~ 2 / D . Again, the proof is based on the similar argument as in the ID case, 
except that one has to replace the characteristic length Iq by the volume Iq in arbitrary 
dimension. The numerical confirmation of the scaling for the 2D case is shown in the 
inset of figure [2} 

For comparison, we plot in the same inset the ground state energy obtained from a 
random-amplitude disorder (i.e., one where to each site corresponds a random potential 
Vi, uniformly distributed in the interval [0, V']). We see that the rate of convergence 
of the energy for a Bernoulli distribution is quicker than for a uniform distribution. 
Because of this, the potential with a Bernoulli distribution is ideal for ground state 
energy convergence as well as localization of low energy states. This may be quickly 



seen as follows, while further details are given in Appendix A If the potential takes 
value zero with a positive probability p, one can define islands of zero potential as the 
natural locus of the ground eigenstate. In cases where the values of the potential are 
positive, even when they are arbitrarily close to zero, the low potential islands can still 
be defined, using a (volume-dependent) energy cut-off, but the contribution to energy 
from such islands will in this case be larger. Accordingly, the convergence of the ground 
state energy to zero will be faster in the former case. Among the distributions which 
assign a fixed probability p to the zero value of the potential, the optimal ones to work 
with are the Bernoulli distributions, in which some positive value V is assumed with 
probability 1 — p. While the rate of convergence of the ground state energy to zero 
is comparable for all such distributions, the advantage of the Bernoulli distribution is 
that it localizes the low energy states more clearly. The analysis of other potential 
distributions is more complicated, because it necessitates introducing an energy cut- 
off to define low potential islands. The Bernoulli potential is also easier to work with 
numerically, since it requires sampling fewer potential realizations. 

Although we do not expect any of the properties discussed in the following to depend 
on the specific type of bounded disorder, we clearly see that the Bernoulli choice yields a 
much faster convergence of the ground state energy to the desired asymptotic behaviour. 

3.2. Lif shits tail 

Since the disordered potential is chosen to be finite, the spectrum of the system is 
bounded from below. As a consequence, the spectrum is expected to exhibit a so- 



called Lifshits tail [10 25 , in the sense that the cumulative density of states (cDOS) 
jV(e) = j e de' V(e') behaves at low energies as 

M{e) ~ exp [-c(e - V^)~ D ' 2 ] . (3) 
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Figure 2. Properties of non-interacting dirty bosons in 2D. Main figure: the 
cumulative density of states Af(e) shows a Lifshits tail (purple dashed line) given by Q . 
Here we have taken a 200x200 lattice, and averaged over 10 realizations. Inset: ground 
state energy Eq in 2D plotted versus linear size of the system L; results for Bernoulli 
disorder are plotted in blue and scale as ~ 1/ log L, while results for random-amplitude 
disorder are plotted in orange and scale as ~ log log LJ log L. 



The cDOS and the associated Lifshits tail for the considered 2D system were obtained 
numerically and are shown in figure [2} The Lifshits states lie on well-separated islands, 
and have almost degenerate energies since the islands' diameters are approximately the 
same. 



4. Interacting case 

In this Section we will present two common theoretical approaches used to describe a 
Bose gas with short-range interaction potential. The first one corresponds to a simplified 
version of Multi-Orbital Hartree-Fock (sMOHF) method, based on an expansion into 
single-particle eigenstates. As we will see, sMOHF is suitable to describe a weakly- 
interacting system, whose ground state occupies few disconnected large islands of zero 
potential. For more sophisticated version of the Hartree-Fock approach applied to 



bosons, see 21 . The second method is the standard Gross-Pitaevskii (GP) equation, 
which describes the dynamics of a fully-coherent matter wave, and correctly describes 
also the regime of strong interactions, where large overlaps between the islands and 
self-interaction on each island play an important role, overlap between the islands. 

Before going into details, let us estimate when the sMOHF description based on 
non-interacting single particle states should be valid. As we have discussed, this regime 
corresponds to the Lifshits glass region. Let us assume that we have r filled islands of 
the characteristic size Iq, and "volume" Iff, so that 

riff = CL D , 

where the proportionality constant C depends, in general, weakly on p and the density 
p, since it results from the complex interplay between the kinetic and interaction energy. 
In the following we will ignore this dependence. 
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The Lifshits glass regime occurs then for an inter-particle interaction coupling g 
smaller than the characteristic value g^, at which the kinetic energy is comparable to 
the interaction energy. An estimate of the characteristic coupling is easily obtained by 
equating the kinetic energy per particle ~ with the interaction energy ~ gN/fl^ . 
Using the above expression for f , we obtain 

g ch ~ C/pll ~ C(| logp|) 2 / D /p(| \ogL\f' D . (4) 

This scaling provides a good estimate of the coupling at which the two energies depicted 
in figure [3] start to diverge. 

4-1. Multi-orbital Hartree-Fock approach 

In the sMOHF treatment, the wave function is expressed in the single particle 
eigenstates, and is taken to be in the product form 

v ^ 

where |0) is the vacuum of the system and a\ creates a particle in the kth non- interacting 
eigenstate (i.e., orbital). We consider repulsive interactions of strength g > 0, which 
yield an interaction energy given by 

(Uint) = I f dx ^(x)^(x)^(x)V>(x) = 



2 _ 

2 / dx S^k 1 ( x )^k 2 (x)^k3(x)^k4(x)4 i 4 2 a k3 a k4 , (6) 

where the sum runs only over terms conserving number of particles, i.e., such that 
ki + k 2 = ks + ktt. The diagonal terms (alala k a k ) contribute a factor | nk(nk — l)Ok,k, 
where we have defined the overlaps Ok,j = f dx|-^ fe (x)| 2 | , :; (x)| 2 . 

There are two other possibilities which conserve the number of particles: (k^ = k\ 
and kj, = k 2 ), or (fc 4 = k 2 and k% = ki). These two terms, called direct and exchange 
terms, are equal for bosons, and their sum contributes a factor | J2j^k^ n k n jOk,j- 

The average interaction energy reads then: 



(7) 



nk{n k - 

The occupation probabilities n k yielding the ground state may now be found by 
minimizing the total energy 

E = J2 n kE k + (U int ), (8) 

k 

subject to the constraints of normalization and positivity of all n k , 

J2n k = N, Vk,n k >0. (9) 

k 



We present the analytic solution of this problem in Appendix B 
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Figure 3. Ground state energies obtained by sMOHF (dashed) and GP (continuous) 
for various disorder densities p. The value of the coupling constant at which the 
energies obtained by the two methods start to diverge is estimated by equation Q. 
Results obtained by averaging over 5 simulations and V = 5t. 



4-2. Gross-Pitaevskii approach 

The properties of the system may be analysed also in the framework of the usual GP 
equation, 



which describes the behaviour of an interacting Bose gas in terms of a single-particle 
coherent wave-function. We find the ground state and the chemical potential /j, by 
imaginary time evolution, applying an operator split method. The ground state energy 
E is related to fi by the formula 



5. Comparison of sMOHF and GP approaches 

In this Section, we present the numerical solutions of the sMOHF and GP equations, 
comparing the performance of the two methods. In particular, we discuss and compare 
the ground state energy, the superfluid fraction, and the fractal dimension. As we will 
see, in the regime of strong disorder (Lifshits regime) sMOHF provides ground state 
energies in accord with GP. Nonetheless, sMOHF has the advantage of being insensitive 
to convergence issues, and its solution is generally much faster than GP. The fact that 
in this regime both methods agree confirms the intuition that the particles populate 
low-lying eigenstates, the Lifshits islands. In the regime where interactions dominate 
over the disorder, we will show instead that GP provides energies which are perceptibly 
lower than sMOHF. Unless otherwise noted, the numerics presented in this section have 
been obtained for a system on a 2D lattice with 32 x 32 sites, iV part = 10 4 particles, a 
Bernoulli potential with V = 5t or V = 50t, and we have set t = h 2 /2m = ks = 1. 
Where needed, we have performed appropriate averages to obtain results which are 
independent of the particular disorder configuration. 




(10) 




(11) 
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a) 




Figure 4. Superfluid fraction p s { as a function of potential density p and interaction 
strength g for a) V = 5t b) V = 50t. 

5.1. Energy 

In figure [3] we compare the energies obtained through the sMOHF approach with those 
obtained from the GP equation as a function of the interaction strength g. For a coupling 
constant g smaller than a characteristic value g c h, the minimization of expression (j8l) 
and the GP equation yield the same values of energy proving that sMOHF method 
correctly describes the system in this limit. The characteristic value g c ^ is given by the 
interaction strength at which the energy obtained from the sMOHF starts to differs from 
the one obtained by the GP equation. The range of g for which sMOHF approach is an 
equivalent description shrinks with increasing p, as the scatterers become increasingly 
sparser, leaving large regions of zero potential. The range of applicability of sMOHF also 
shrinks with decreasing V. For values of the interaction g ^> g c ^, the disorder plays a 
negligible role and the energy per particle saturates to the analytic value E = E + gp/2, 
where E is a constant that depends on V and p. This analytic result is recovered by 
the numerical solution of the GP equation, but not by the sMOHF approach. 



5. 2. Superfluid fraction 

The superfluid fraction p s f may be defined in terms of the energy change of a periodic 
system in response to twisted boundary conditions along one direction. To calculate 



this quantity we solve the GP equation (10) with boundary conditions ipi+Lj — e r»j> 
and we obtain 

2mL 2 g(g) - E(0) 

P« = -tf $2 ' ( 12 ) 

where E(<&) is the energy per particle of a system with the total phase shift $. Our 
results are depicted in figure |4j In the regime where sMOHF provides a good description, 
the SF fraction is negligible because the localized single-particle eigenstates experience 
an exponentially small energy change due to the imposed phase shift. The GP equation 
instead yields a perceptibly lower energy than sMOHF when the superfluid fraction p 8 f 
becomes sufficiently large (p s f > 0.1). 



Glass to superfluid transition in dirty bosons on a lattice 



11 



lo gi P 

1.0 



lo gi P 
2.00 r -. 



0.7 1.95 
1.90 
06 1.85 



1/logL 



0.1 0.2 0.3 
0.00 0.05 0.10 0.15 



••. o 
0.20 



0.25 



1/logL 



Figure 5. Scaling of the logarithm of the participation number log L P versus the 
inverse of the logarithm of the lattice size. The slope of the lines, being the logarithm 
of c, allows for determining the fraction of space occupied by the state. Main figure: 
ID lattice of length L, and interaction strengths g = O.l(empty) and g = 10(filled). 
The potential densities are respectively p = 0.4(red circles) and 0.7(blue squares). 
Inset: 2D lattice of linear size L, interaction strength g — 10, potential parameter 
p = 0.6(green diamonds), 0.8 (purple triangles). 



5.3. Fractal dimension 

As we have discussed in the previous sections, the ground state of a weakly-interacting 
fluid is localized on one or a few Lifshits islands. When interactions start to play 
an important role, the extension of the ground state increases generally until the gas 
occupies all available space. In order to provide a quantitative measure of the extension 
of the ground state, we calculate its fractal dimension d* [j], and introduce here the 
concept of fractional occupation c. 

The fractal dimension is defined as a minimum d* such that 
P 

lim — = c, c> 0. (13) 

Here P = 1/ f dx|-?/> (x)| 4 = 1/O 0) o is the so-called "participation number" of the 
ground state fij. For an extended state such as a plane wave one finds that P is 
proportional to the volume of the system, and therefore the fractal dimension equals 
the Euclidean dimension, d* = D. For a state which is instead localized in a volume 
Q, the participation number behaves as Q, and the corresponding fractal dimension 
vanishes if Q grows slower than any power of L. For instance, the ground state of the 
non-interacting system is localized on an island of volume Q ~ logL which for all a 
grows slower then L a . So, in general, one has < d* < D. 

Analysing the results obtained from both sMOHF and GP approaches, we find 
that d* = D for any g > 0. This may be seen as follows. Assuming that d* is bounded 
strictly below D, we see that the interaction energy (gN/2) J dx |?/>o(x)| 4 will diverge. 
Physically, there is not enough space for the particles to keep the interaction energy 
bounded unless they spread out through a non-zero proportion of the whole space. If 
log L P = D — f(L), and lim^oo L D / L d * = c, then it must be that L^ L ^ — > c, so 
f(L) ~ log(c)/ log(L). Numerical results confirm this behaviour in both ID and 2D. To 
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Figure 6. Fractional occupation c of the ground state as a function of potential density 
p and interaction strength g for a) V = 5t and b) V = 50i. Here we have considered a 
2D lattice and c was obtained by averaging over 5 potential realizations. 



show this, in figure [5] we plot log L P as a function of 1/ logL for a ID system and a 2D 
system (inset). 

The quantity c, which takes values between and 1, can be interpreted as the 
fraction of space occupied by the ground state, and we will therefore refer to it as 
the fractional occupation of the ground state. We show our results for the fractional 
occupation in figure [6j 

For strong disorder (V 3> t) the fractional occupation c displays at intermediate 
interactions a plateau at ~ p. This coincides with a similar plateau in the plot of the 
superfluid fraction (cf. figures |4jo and^). In ID the explanation for such behaviour of c 
and p s f is the following. For small g the wave function fills only a part of zero-potential 
space. Indeed, for g = 0, we have a linear operator and therefore the ground state 
is approximately the wave function on the longest island of zero potential, yielding a 
nearly zero c. For g small such that gN/2 is approximately l/(logL) 2 , the ground state 
spreads to many of the long islands: the kinetic energy increases a little (on the order 
l/(logL) 3 ) while interaction energy decreases by a factor l/(# of islands). As g grows 
to be on the order of 1, it supports itself on all sites of zero potential save for the shortest 
islands, so c ~ p. This is where the graphs level off in (Fig [6]), where the ground state 
does not have large norm on short islands of zero potential and sites of V potential. 

For large interactions one may use the Thomas- Fermi Approximation, cf. (SI. 



Suppose we choose an ansatz such that the wave function equals on sites of 
zero potential and ^]j^ 2 , then energy optimizing m would be m 2 = p + and 
1 — m 2 = q — If gp < pV, then the above ansatz is invalid, and the ground 

state stays exclusively on sites of zero potential. If we compare the energy of a wave 
function equal to ^=j= on sites of zero potential and zero elsewhere to a wave function 
that is equal to -7= everywhere, and ignore the kinetic energy terms (which we can 
control by another parameter), we see that the first wave function provides a better 
variational ansatz if ^ < pV. This means that for gp < 2pV, we are in an area 
where c p, that is, the kinetic energy and interaction energy are balancing, but any 
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significant support of the ground state on the sites occupied by scatterers is too costly. 
Since the wave function does not spread significantly in this regime, also the superfluid 
fraction remains constant. 



6. Comparison of temperature and interaction effects 

In this section we turn to the analysis of the effects of non-zero temperature in a non- 
interacting system. We will see that those remind very much effects of interactions at 
zero temperature, and yet are quite different. 



6.1. Single-particle density distribution 

For free boson systems, the effects of temperature may be studied by expanding arbitrary 
states in the basis of non-interacting eigenstates, in analogy with what we did to treat 
the interacting case. In this approximation, the average density of the gas reads 

p(x) = L- D 5> fc |<Mx)| 2 . (14) 
An approximation of a typical wave-function is then given by 

^(x) = ^v^e^*V*(x), (15) 

where e l ^ k are random phases, and w k is a Gaussian random variable with distribution 
~P(wk) = n^ 1 e~ Wk ^ Uk , such that (w k ) v = n k . 

The occupation factors depend on the statistics of the particles. Identical bosons 
follow the Bose-Einstein distribution 

Uk = e(^J/T-r ( 16 ) 

Distinguishable (or classical) particles would instead populate the non-interacting 
eigenstates following a Boltzmann distribution, n k = e~( Ek ~^' T . 

At zero temperature, only the lowest energy state will have a non-zero population, 
i.e., n k = N5 0t k- Non-zero temperatures will modify the occupation probabilities, 
redistributing population over the higher energy states, but we expect that for T <C t 
only the lowest energy states (the Lifshits states) will be populated. 

We note here that in an interacting 2D Bose gas of finite extension one 
generally expects a crossover between two qualitatively different regimes. At very low 
temperatures (T < Tq) phase correlations in the gas decay algebraically with distance. 
For temperatures larger than the critical value Tc instead, phase correlations decay 
exponentially with distance. This switch in the decay of correlations can be traced back 
to the dissociation of bound vortex-antivortex pairs at T > Tc- The crossover becomes 
a real phase transition in an infinite system, and the underlying mechanism goes under 



the name of Berezinskii-Kosterlitz-Thouless (BKT) transition 26 . 

Our approach to finite temperatures is rather crude, and is not well suited to 
describe BKT physics, which requires considering both temperature and interactions at 
the same time. However, for the range of parameters considered in this manuscript, our 
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Figure 7. Comparison of temperature and interaction effects. First two rows: 
thermally averaged single particle density with respectively Boltzmann (top) and Bose 
(centre) statistics; from left to right, T/t —0, 0.1, 0.2, 0.5. 

Bottom row: single particle density of an interacting gas; from left to right, gp/t = 
0, 0.1, 2, 10. 

The red dots indicate the sites occupied by the disordered scatterers. The average 
density is p — 0.9particles/site. 

model provides a qualitative, although simplified, explanation of the complementary 
roles payed by temperature and interactions. Moreover, a similar scenario to the 
one considered here is present in three dimensions, where a true BEC exists at low 
temperatures. 

6.2. Non-zero temperature 

As temperature is raised, the cloud gets to occupy more and more islands; when T < V 
it will populate all unoccupied sites, and for T 3> V it will occupy even sites inhabited 
by scatterers. 

The averaged densities in presence of positive non-zero temperature are shown in 
figure [7^) and b). While a classical cloud quickly spreads over various low energy states, 
a bosonic gas at low temperatures tends to condense in a single or few Lifshits states. 
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120 



Figure 8. Occupations of the single particle states for the simulations shown in 
columns 2,3, and 4 of figure [7J The lines are: Boltzmann (yellow dotted), Bose (red 
dashed), and sMOHF (blue continuous). Results averaged over 200 simulations, x-axis: 
index labeling the single-particle state; y-axis: occupation. 



Interaction effects in a bosonic cloud at T = 0, instead, yield at first sight a state which 
very much resembles the Boltzmann case. Even for moderate repulsive interactions, the 
cloud spreads over various Lifshits states. Effects of interactions at T = are shown in 
figure [7}:. 

The occupation probabilities nk for the three cases analysed above (Boltzmann, 
Bose, and interacting) are shown in figure [8j Since the lowest-energy single-particle 
states correspond to the localized and well separated Lifshits states, a population 
distribution which is narrow in energy corresponds to a ground state localized on one 
or very few well-separated islands, the Lifshits glass. As the temperature increases the 
Bose liquid tends to have a rather narrow population distribution, while the Boltzmann 
and interacting distributions quickly spread. Nonetheless, a few important differences 
may be noted comparing the two latter maybe seen in Figs. [7£i and ^jp. While 

the population distribution for the Boltzmann case smoothly decreases, with a long tail, 
the distribution for the interacting case has a local minimum at very low k, while at 
larger k goes through a maximum and then quickly drops to zero. These features can 
be understood in the following way. The eigenstates with lowest energies tend to have 
higher inverse participations (overlaps), i.e., the Lifshits states are rather concentrated 
on a single or few islands while states at intermediate energies are delocalized over the 
whole system. For an interacting system is therefore preferable to occupy states indexed 
by states with intermediate k values, as their larger spatial extension helps to reduce 
the interaction energy. At very high energies (k ^> 1) instead the states have again very 
high Okk, since they become completely localized on single sites, and their populations 
abruptly drop to zero, as it is too costly energetically to put particles with repulsive 
interactions in small regions of space. 

7. Conclusions 

In this paper we have considered a bosonic lattice gas in the presence of Bernoulli 
disorder, given by randomly-localized identical scatterers. We have shown that in this 
case it is possible to provide precise analytical estimates even for the interacting case. 
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We have compared two theoretical schemes, the simplified multi-orbital Hartree-Fock 
and the Gross-Pitaevskii approaches, showing how the first is very accurate in the 
glassy regime of strong disorder, but it fails when interactions bring the system into an 
extended state and the superfluid fraction reaches values p s f > 0.1. 

Further, we have shown analytically that the fractal dimension for this kind of 
potential tends to the physical dimension of the system. As a result, we have introduced 
a quantity termed fractional occupation, which characterizes the typical extension of the 
system. When the latter becomes of order p, i.e., the fraction of physical space where 
scatterers are absent, the system crosses over from the Lifshits to the Bose glass. Finally, 
we have discussed the similarities and differences between effects due to interaction and 
temperature. This question is of fundamental importance for ongoing experiments, since 
in a laboratory the two effects are often difficult to distinguish. We hope that our results 
will provide new insights in the complex route towards understanding the interplay 
between disorder, interactions, and temperature in quantum mechanical systems. 
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Appendix A. Comparison of the Bernoulli and uniform distributions 

The choice of potential may change the behaviour and scaling of the ground state energy 
in a single-particle system. We will consider potentials that are bounded below at and 
to simplify notation let P(x) = P[Vj < x] for a given distribution. Any potential that is 
bounded below can be shifted so that the bottom of its support is zero. Then P(x) = 
for x < 0. For the ground state, we are interested in the distribution around zero; 
the behaviour of the distribution outside of the neighborhood of zero does not affect 
limiting behaviour. Assume that the probability of the potential being less than e in a 
neighborhood of zero is given by P[V < e] = a + me 7 , where a G [0, 1), m > 0, and 
7 > 0. Then the expected potential energy of a ground state supported on the island 
with each site having potential less than e is: 
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m7£7+1 (A.l) 



2(7 + l)(a + me7)' 

For a specific e and related probability P(e), the typical largest island composed 
by sites where the potential is less than e has size 
^ log(P(e)(l-P(e))L) 

Lmax te^pje) (AJ) 

Then the energy of the ground state is expected to be approximately 
(|logP( e )|) 2/D | mje^ 



(log(P(e)(l-P(e))L)) 2 / D 2( 7 + l)(a + me 7 ) 

For this to be minimized in the large system limit, the order of each term must be 
approximately the same, so e must scale with L as: 

^ * (i^F 5 ' (A ' 4) 

Plugging this into the kinetic term should give an approximate rate of convergence 



(A.5) 



(\ogL) 2 / D 

For the Bernoulli distribution, a = p and m = 0, so the rate of logL~ 2 / D is 
recovered. For the uniform distribution, a = and m = yj and has rate of convergence 
( ^lo'gf L ) 2//D which is slower than the Bernoulli distribution. In fact, any distribution with 
a 7^ 0, which means the potential can equal exactly zero with positive probability, must 
have its associated ground state energy converge to zero faster than for distributions 
with a = 0. Considering distributions with a > and comparing the cases m = and 
m > 0, the ground state energy in the former case will be larger than in the latter case, 
but asymptotically they will converge to zero at the same rate. More importantly, in 
the case where m — 0, the longest island is clearly defined whereas its definition involves 
an arbitrary (or artificial) cutoff, in the case where m > 0. 

Appendix B. Solution of the sMOHF equations 

The occupation probabilities n k may be found by minimizing the total energy 



k 



(B.l) 



where O = 10 — Diag(O), subject to the constraints: 

J2n k = N, Vk,n k >0. (B.2) 
k 

The minimization yields the linear system (one equation for each k): 

\i + Ik - E k O kyk 

n j°j,k = + -s-, (B.3) 
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where Lagrange multipliers 7fc and /x correspond to the conditions of positivity of each 
rife and global normalization, respectively. We introduce the vector 



v k 



I 1 — E k Ok,k 

g ~2~ 



which allows us to write the solution as 
n 



(B.4) 



(B.5) 



9 

Now, 7 and /x should be chosen such that the following conditions, named KKT [27] , 
are fulfilled: 



7 fc > & 7 fc n fc = & 2J n k = N. 



(B.6) 



The KKT conditions allow for solving the optimization problem in the case in which 
some of the constraints are in the form of inequalities. These conditions lead to a 
nonlinear system for 7 and /x. 
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